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ABSTRACT 

Magnetically confined mountains on accreting neutron stars are promising sources 
of continuous-wave gravitational radiation and are currently the targets of directed 
searches with long-baseline detectors like the Laser Interferometer Gravitational Wave 
Observatory (LIGO). In this paper, previous ideal-magnetohydrodynamic models of 
isothermal mountains are generalized to a range of physically motivated, adiabatic 
equations of state. It is found that the mass ellipticity e drops substantially, from 
e « 3 x 10 -4 (isothermal) to e « 9 x 10 -7 (non-relativistic degenerate neutrons), 
6 x 10 -8 (relativistic degenerate electrons) and 1 x 10 -8 (non-relativistic degenerate 
electrons) (assuming a magnetic field of 10 12,5 G at birth). The characteristic mass 
M c at which the magnetic dipole moment halves from its initial value is also modified, 
from M c /M « 5 x 10" 4 (isothermal) to M c /M « 2 x 10" 6 , 1 x 10" 7 , and 3 x 10" 8 
for the above three equations of state, respectively. Similar results are obtained for 
a realistic, piecewise-poly tropic nuclear equation of state. The adiabatic models are 
consistent with current LIGO upper limits, unlike the isothermal models. Updated 
estimates of gravitational- wave detectability are made. Monte Carlo simulations of the 
spin distribution of accreting millisecond pulsars including gravitational- wave stalling 
agree better with observations for certain adiabatic equations of state, implying that 
X-ray spin measurements can probe the equation of state when coupled with magnetic 
mountain models. 
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1 INTRODUCTION 

Neutron star spins in low-mass X-ray binaries (LMXBs), 
measured from X-ray pulsations or thermonuclear burst 
oscillations, a r e fou n d to lie in the r ange 95 — 619 Hz 
(IChakrabartvl |2008l: iGallowavl 120081 : IWatts et al.1 120081 : 
Gall oway et al.1 l2010h . The upper end of this range falls 
well short of the c entrifugal breakup frequency for most 
equations of state (|Cook et al.l 1 19941 : lHaensel et al.l [l999; 
IChakrabartvi r2008). even though the objects accrete enough 
angula r momentum during their X-ray lifetime of 10 7 — 10 9 
years (IPodsiadlowski et al. l 20021) to spin up to 1.5 kHz < 
Us < 3 kHz (|Bildstenlll998l : lchakrabartv et alJl2003h . This 
discrepancy cannot be attributed to an observational se- 
lection effect, because the Rossi X-ray Tim ing Explorer 
(RX T E) remains sen sitive up to 2 kHz ( Chakr abartvl 
120081 : IGallowavl l2008h . To desc r ibe t he apparent spin 
clustering and cut-off, iBildstenl (|l998l) invoked gravita- 
tional radiation torques to stall the spin-up process; see 
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also IPaoaloizou fe Pringlel (| 19781 ) and IWagoneH (| 19841 ). To 
achieve this, a mass quadrupole moment of order ^ 2 x 
10 38 g cm 2 is required. 



Quadrupoles on 
two kinds: (i) core 



accreting neutron stars are of 
deformations, e.g. from r-modes 



(iBrink et al.ll2004l ; iNavvar fe OweJ200d Bondarescu et al.1 
l2007l ) and (ii) per manent crust al deformations, e.g. sup - 
ported by the rmal (IBildstenl 1 1 998l: lUshomirskv et al.ll2 000) 
or m a gnetic (|Brown &; Bildsten| Il998l: iMelatos &; Phinneyl 



20011: IChoudhuri fc KonaJ l2002inPavne Melatod 1200 ^ 
IVigelius"fc Melatos] 120081 ) gradients. In the absence of 
a magnetic field, the maximum crust al quadrupole de- 
pends on the bre aking strain ([Ushomirskv et al ] l200d : 
lHaskell et all l2006h and can be as large as ~ 10 g cm 
in the light of recent molecular dynamics simulations 
([Horowitz &: Kadaul I2OO9I ) . When magnetic stresses are in- 
cluded, the quadrupole increases, as matter is funnelled to 
the magnetic pole s of the star and co mpresses the magnetic 
field l aterally (lHameurv et al. _ 19831; [M elatos Phinne^ 



20011: IChoudhuri fc Konarl 120021 : IPavne Melatod 120041 : 
Vigelius fc Melatod 120081 ) 
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IPavne fe Melatos! (|2004l ). hereafter PM04, calculated 
self- consistent, axisymmetric, ideal- magnet ohydro dynamic 
(ideal-MHD) equilibria of isothermal magnetic mountains as 
a function of accreted mass M a . They found that the mag- 
netic field distorts appreciably for M a ^ M c ~ 1O _5 M0, 
in accord with the ph enomenological field decay relation of 
IShibazaki et al.l |l989) and well above previous calculations, 
which predicted M c ~ 10 _10 Mq without including the 
back reaction from the compressed equatorial magnetic field 
(lHameurv et al.lll983l: iBrown fc Bildstenl[l998l :l Lit win et al.l 
2001). Payne & Melatos (2007) showed that the mountain 
oscillates stably in a superposition of Alfven and acous- 
tic modes when perturbed, following a transient adjust- 
ment via t he undular submode of the magnetic buoyancy 
instability (iMouschoviad 1 19741: iHughes fc Cattaneol Il987l : 
IVigelius fc Melatod I2008T ). I Vigelius fc Melatod (|2o"o8h found 
that the equilibrium state remains mountain-like after this 
transient instability, with the mass quadrupole moment de- 
creasing by ~ 30 per cent. Ohmic dissipation contributes 
to the decay of the mass quadrupole by allowing slippage 
of accreted matter across magnetic field lines, with a re- 
sistive relaxation timescale of 10 5 — 10 8 yr d epending on 
the co nductivity (|Vigelius fe Melatod l2009bh . IWette et al.1 
(2010) examined the subsidence of mountains into a fluid 
crust, generalizing earlier calculations on a rigid surface, and 
found that the quadrupole shrinks by up to « 60 per cent. 

The existing literature on magnetic mountains, sum- 
marized above, suffers from several limitations. First, the 
time-dependent feedback bet ween the magneto-sphere and 
the accretion disc is neglected jRomanova et al. 2003, 2004; 
iKulkarni fe Romanoval I2QQ8I : lLong et al.1 120081 ). Secondly, 
the mountain should solidify into a body-centred-cubic crys- 
tal as it sinks, when the ionic coupling parameter exceeds 
the c rystalliz ation thre shold (iFarouki fe Hamaguchil Il993l : 
iHorowitz k, Berrvl 120091 ). This occurs at different depths, 
depen ding on the l ocal composition, density and temper- 
ature (jBrownl 2000). The sudden transition to a solid af- 
fects the magnetic line-tying boundary condition, which 
now depends on the local magnetic stresses and criti- 
cal strain. Thirdly, a nuclear reaction network that fol- 
lows accreted matter el ements as they descend ha s not 
yet been implemented (lHaensel fe Zdunikl Il990al lbl. 120031 : 
Cha mel fe Haensell l2008h . Deep crustal heating deposits 
1.5 — 1.9 MeV per accreted baryon (lHaensel k, Zdunikl 
2008), reduces the Ohmic decay time-scale, and intro- 
duces thermal and electri cal conductivity gradients due to 
compositional variations (|Chamel fe Haensell l2Q08h , all of 
which affect the mountain structure. Finally, the equa- 
tion of state (EOS) of the accreted matter needs to be 
modelled realistically. The calculations cited in the pre- 
vious paragraph all utilize an isothermal EOS, an accu- 
rate model for very low mas s mountains with maximum 
density p max ^ 10 6 g cm" 3 (jShaoiro fe Teukolskvl [l983h . 
The isothermal EOS is too soft and does not accurately 
represent all pressure components (e.g. degenerate neutron 
and electron pressures in the inner crust) for realistically 
sized mountains with p ma x 5. 10 14 g cm~ 3 °r equivalently 
M a < 10~ 2 M^ (iShapiro fc Teukolskvl Il983l : iBrownl l200ol : 
IChamel fe Haensell l2008l ~ 

This work aims to quantify how the EOS influences the 
structure of the magnetic mountain and its mass quadrupole 
moment. It turns out that the effect is large. In Section [21 



we generalize the Grad-Shafranov framework for solving nu- 
merically the MHD equilibrium problem to incorporate an 
adiabatic EOS. The numerical algorithm is validated against 
published isothermal results in Section [3] We directly com- 
pare the structure of adiabatic and isothermal magnetic 
mountains in Section [4] quantifying the relation between 
the accreted mass and measurable quantities such as dipole 
moment and ellipticity. In Section [5] we approximate the 
realistic EOS in the neutron star crust by an effective poly- 
trope and calculate the structure of the associated mountain. 
In Section [6j we examine the implications of the theoretical 
models for gravitational- wave (GW) stalling of LMXB spins. 
The detectability of magnetic mountains as GW sources is 
assessed briefly in Section [71 revising the latest estimates in 
Vigelius fe Melatos (20093). 



2 HYDROMAGNETIC EQUILIBRIUM 

To compute the structure of a magnetic mountain with 
an adiabatic EOS, we generalize the isothermal Grad- 
Shafranov solver described in PM04 to handle a general, 
barotropic, pressure-density relation of the form P(p) = 



Kp^ 1 = Ap^ , wUer e n is trie polytropic index and 1 is 
the adiabatic index (|Paczvnskil[l983l : IShapiro fe Teukolskvl 
19831 ). 



l + l/n 



2.1 Grad-Shafranov equation 

Let us define a spherical coordinate system (r, 0,0), where 
6 — is the magnetic symmetry axis before accretion begins 
and the neutron star surface is situated at r = Ri n (i.e. the 
inner boundary of the simulation; see Appendix [A} . Time- 
dependent ideal MHD and resistive simulations of magnetic 
mountains in ZEUS-MP show that the magnetic field relaxes 
to an almost axisymmetric configuration (deviation from ax- 
isymmetry < 1 per cent) within a few Alfven times, follow- 
ing a transient, Parker- type instability ( Vige iius &; Melatos! 
2008). Hence, to a good approximation, the magnetic field 
is given everywhere by 



x e^, 



(i) 



where ip(r, 0) is a flux function. In the steady state, the MHD 
equations reduce to 



VP + + (A 2 V0 VV> = 0, 

where A 2 denotes the Grad-Shafranov operator, 



47rr 2 sin 



+ ■ 



sin# d 



1 d 



2 80 \ sinOdO 



(2) 



(3) 



We solve the projection of equation ([2]) along B by the 
method of characteristics. The result depends critically on 
the EOS. Under isothermal conditions, i.e. P — c^p, we find 



dF(VQ 
d^ 



exp[-(0 - o )/cs], 



(4) 



where 0o denotes the reference gravitational potential at 
the neutron star surface, and c 2 is the isothermal sound 
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speed (jPavne fe Melatos 20041 ). Under adiabatic conditions, 



i.e. P = Kp , we find 



jgw L _ (r- i)(0-0o 



d^ | r^/r^^jcr-i)/! 



i/(r-i) 



(5) 



The pressure along a flux surface -0 under isothermal and 
adiabatic conditions is given by 

P = FWexp[-(^-^)/c s 2 ], (6) 



and 



(r-i)(0-0p) 

r^!/r[_F(^)](r-i)/r 



r/(r-i) 



(7) 



respectively. Formally speaking, F(ip) is an arbitrary func- 
tion of the magnetic flux in equations (H])-©. Equation 
(|6]) is the usual barometric formula; the base pressure F(ip) 
varies from field line to field line, and P decreases with arc 
length along any particular field line because \(f)\ is inversely 
proportional to r. Equation (0 behaves similarly, but its 
form is not barometric, in the sense that F($) does not fac- 
torize out. 

In order to establish a one-to-one mapping between the 
initial (pre-accretion) and final (post-accretion) states that 
preserves the flux freezing encoded in the mass-continuity 
and magnetic-induction equations of ideal MHD, we require 
that the final, steady-state, mass- flux distribution dM/dip, 
defined as the mass enclosed by the infinitesimally separated 
flux surfaces ip and ip + dip, equals that of the initial state 
plus the accreted mass. This approach uniquely determines 
F(ip) through 



W ' 2tt dip 



L 



ds r sin#| V^| 1 e 



-{4>-4>o)/ci 



for the isothermal EOS and 

/ \ r 
K I r\i\/r y 



AM 

Aip 



(8) 



(9) 



/ ds rsin^lV^I 1 
Jc 

(r-i)(0-0 o ) 



i/(r-i)' 



for the adiabatic EOS. This approach is self-consistent and 
therefore preferabl e to guessing F(ib) (lHameurv et al.l ll983: 
iBrown fe Bildstenl Il998l : iMelatos fe Phinnevl l200lh but it 
renders the solution more difficult. [Duez Mathisl (|2010h 
also solved self-consistently for F(ip) by minimizing the to- 
tal energy while conserving invariants like the helicity and 
mass-flux ratio.] The integrals in equations J5]) and (0 are 
performed along the magnetic field line ip = constant. In ac- 
cordance with earlier work, we prescribe the mass- flux dis- 
tribution in one hemisphere to be 



M(tp) 



M a [l -exp(-V>M0] 



(10) 



2[1 - exp(-b)] 

where M a is the accreted mass, -0* labels the flux surface 
emerging from the magnetic equator, ip a labels the field line 
that closes just inside the inner edge of the accretion disc, 
and we write b = ip*/ip a . Equation ([TO]) ensures that « 63 



per cent of the accreted mass accumulates within the polar 
cap ^ -0 < ip a for < V>* • 

The gravitational acceleration is assumed to be con- 
stant in this paper, with a gravitational potential of the 
form <j)(r) — GM*r / Rf n . This assumption is justified, be- 
cause the mountain never rises more than ~ 10 4 cm above 
the hard surface at r = R in (see Section [43)) . A simple nu- 
merical check shows that the altitude above r = R in where 
the density distribution falls to zero changes by ~ 2 per cent 
when <j)(r) = GM*r/Rf n is replaced by 0(r) = -GM*/r. 
Self- gravity is also ignored, although the correction M a /M* 
to the gravitational potential is significant in LMXBs with 
Ma>lO _1 M . 

We conduct our numerical simulations as follows: a fixed 
dipolar magnetic field at the inner radial boundary of the 
numerical mesh is assumed, and a prescribed amount of ac- 
creted matter M a (described by one of the EOS in Table 
[J) is added into the simulation volume according to the 
mass- flux relation (|10p . We then allow the system to re- 
lax quasi-statically to hydromagnetic equilibrium by solv- 
ing equation Q or ([5]) simultaneously with equation (|5J) 
or (9]) for ip(r,6), using an iterative under-relaxation algo- 
rithm combined with a finite- difference Poisson solver. The 
details can be found in Appendix [A) We adopt the follow- 
ing boundary conditions, as in previous papers (e.g. PM04): 
ip(Rin,0) — ^* sin 2 (surface dipole; magnetic line tying), 
d^/dr(R m , 0) = (outflow), tp(r, 0) = (straight polar field 
line) and dip/d6(r, 7r/2) = (north-south symmetry), where 
Rin ^ r < R m and ^ < tv/2 delimit the computational 
volume. The outer radius R m is chosen large enough to en- 
compass most of the screening currents (isothermal EOS) or 
the outer edge of the accreted matter (adiabatic EOS). 

2.2 Inner boundary 

The nature of the rigid inner boundary at Ri n deserves spe- 
cial mention. It is not the stellar surface; it is not meaningful 
to build a mountain 100 m high and reaching neutron drip 
density at its base on top of a low-density ocean, using a 
realistic EOS. Instead, the outer layers of the neutron star 
are 'constructed' from the accreted material of mass M a . 
Thus, Ri n does not correspond to the neutron star surface 
i?*, but to the depth in the neutron star crust above which 
lies the mass M a (for a given EOS). Since M* and R in are 
fixed, the total mass and radius vary slightly (few per cent) 
between models with different M a but the same EOS (Ta- 
ble [1]) . The inner boundary of our simulation volume Ri n 
represents a solid surface at the corresponding base den- 
sity. This simplification assumes that movement of matter 
below this depth is approximately radial due to compres- 
sion and that the solid-surface prescription is valid. In re- 
ality, accretin g matter is expected to disp lace both radially 
and laterally ([Choudhuri &; Konarl 12002). The lateral flow 
would alter our computed results by decreasing the mass 
quadrupole moment slightly and increasing the magnetic 
dipole moment. We can eliminate this approximation by in- 
jecting t he accreted matter according to the approach advo- 
cated bv llWette et al.l (|2010|), generalizing the latter paper to 
a realistic EOS. Such a procedure is feasible but technically 
difficult; we defer it to fu ture work. 

Referring to fig. 12 of lWette et all (|2010h . the ellipticity 
of an isothermal mountain in the fluid-surface model appears 
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to converge to the saturation ellipticity of the hard-surface 
model as M a increases; the difference in ellipticities relative 
to the hard-surface model decreases from ~ 60 to ~ 25 per 
cent as M a increases from ~ 10 -3 to ~ 10 _1 Mq. We ex- 
pect similar convergent behaviour for adiabatic mountains 
at significantly lower M a , since saturation ellipticities of adi- 
abatic mountains are attained at accreted masses 2-4 orders 
of magnitude below that of the isothermal one (see Fig. 2}. 
Realistic accreted masses in LMXB systems of M a ~ 10 _1 
Mq are 2-6 orders of magnitude greater (depending on the 
EOS) than the accreted masses which we can reliably simu- 
late. At realistic accreted masses, we expect the saturation 
ellipticities of mountains with and without sinking to ap- 
proximately converge. The population-synthesis results in 
Section [6] and GW-detectability estimates in Section [7] de- 
pend solely on the saturation ellipticity. 

2.3 Adiabatic index 

The realistic EOS of a neutron star crust is piecewise adia- 
batic, as discussed in Section [5] However, before modelling 
the realistic EOS, we conduct numerical experiments in Sec- 
tions [3] and [4] to see how the mountain structure depends 
on the adiabatic index Y. In these numerical experiments, 
we employ a purely adiabatic EOS with unique K and Y. 
The values of K and Y are chosen to correspond to den- 
sity regimes of interest in the crust, e.g. degenerate non- 
relativistic electron gas 10 5 < p/(g cm" ~ 3 ) < 10 7 , de gener- 
ate relativistic electron gas 10 7 < p/(g cm -3 ) < 10 12 and 
degenerate neutron gas 10 12 < p/(g cm -3 ) < 10 16 . In an 
ideal electron gas, which is approximately isothermal, ra- 
diation and lattice pressures dominate, but this occurs at 
much lower densities p < 10 4 g cm -3 , which are irrelevant 
to the mountain problem. 

Table [T] displays the magnetic mountain models we com- 
pute here, with the details of their respective EOS. K is 
a function of mean molecular weight per electron, \i e — 
mb/(m u Y e ), according to the scaling K oc /i^ 4 ^ 3 , where vrtb 
is the mean baryon rest mass, m u is the atomic mass unit, 
and Y e is the mean number of electrons per baryon. Under 
the assumption of symmetric nuclear matter, we take p e = 2 
and nib = rn u and hence K is a constant (i.e. independent 
of p). This form of the EOS describes well a completel y 
degenerate, ideal Fermi gas ([Shapiro k, Teukolskvl fl983). 
Hence we use it to model degenerate relativistic electrons 
(n = 3, r = 4/3, K = 4.93 x 10 14 dyn g~ 4/3 cm 2 ), degener- 
ate non-relativistic electrons (n = 3/2, Y = 5/3, K — 3.16 x 
10 12 dyn g- 5/3 cm ) and degenerate non-relativistic neu- 
trons (n = 3/2, T = 5/3, K = 5.38 x 10 9 dyn g~ 5/3 cm 3 ). 



3 VALIDATION IN THE ISOTHERMAL LIMIT 

We assume the following neutron star parameters through- 
out this paper, except where stipulated otherwise: M* = 
1.4M , R in = 10 6 cm, and ip* = 1.6 x 10 24 G cm 2 (with 
-0* = B*Ri n /2, where B* is the polar magnetic field strength 
before accretion begins). The fiducial value of the magnetic 



field, B* 



10 1 



G, is chosen in accord with popula 



tion synthesis mo dels, which predict natal magnetic fields 
of 10 12 - 10 13 G (|Hartman etallll997l ; lArzoumanian et al.1 
2002; iFaucher -Gigu ere &; Kaspi 2006). 



The adiabatic Grad-Shafranov formalism in Section [21 
and the numerical solver described in Appendix [Aj must 
reproduce the results of PM04 in the isothermal limit (i.e. 
n — >• oo, r — >• 1, iT — > c 2 ). In this limit, the isodensity 
contours and magnetic field lines of an adiabatic mountain 
with r — >• 1 must converge to those plotted in figs 4, 5 
and 9 in PM04 for identical accreted masses. As there is 
no unique way to continuously transform an adiabatic EOS 
into an isothermal EOS, we test the adiabatic/isothermal 
correspondence by taking the limit (K, Y) — >> (c 2 , 1) in three 
different ways below. 

(i) We set K — c 2 and let Y tend to unity, such that 
P(p) = c s V- (11) 

(ii) Exploiting the tendency for the surface pressure 
i^surf (0) and density p SU rf (0) to be roughly EOS-independent 
for r « 1, we write P SU rf ~ c 2 p sur f ~ ifpLrf to eliminate K, 
and hence obtain 



p{p) 



pi-r r 

^surf P 5 



(12) 



where P S urf is a function of M a . 

(iii) We take K oc Y and interpolate between a selected 
polytrope (Kq, Yq) and the isothermal target according to 



P(P) 



(r - i)k + (r 



r)c 2 r 
p , 



(13) 



with r 



We apply these three approaches to the case M a = 
1.0 x 10 -5 Mq, starting from a relativistic degenerate elec- 
tron EOS (model C in Tabled]). This EOS prevails over a 
large logarithmic range of densities in a realistic stellar crust 
[10 7 < pi (g cm -3 ) < 10 12 ; see Section[5] and gives way to an 
isothermal EOS in the upper atmosphere (p < 10 4 g cm -3 ). 
We find that all three approaches converge correctly to the 
T = 1 results of PM04 after - 3 x 10 3 iterations. Fig. □ 
displays the mass ellipticity, magnetic dipole moment, and 
grid-averaged residual (relative to the r = 1 result) as a 
function of Y for approaches (i) (red diamonds), (ii) (green 
rectangles) and (iii) (blue triangles). As indicated by Fig. [1] 
the rate of convergence towards the isothermal results dif- 
fers between models. The abnormally high dipole moment 
for case (i) at Y = 1.06 in Fig. [T] is caused by insufficient 
resolution in and can be prevented by scaling the grid log- 
arithmically in 6 to handle the steep magnetic field gradients 
at the equator. We defer this project to future work. 

We compute the mass enclosed within the computa- 
tional grid as a function of iteration number, to track the 
mass lost through the outer boundary. In every converged 
equilibrium, the total mass in the final state is always within 
4 per cent (and typically within 1 per cent) of the initial 
mass. The iterative solver also preserves the divergence- free 
nature of the magnetic field, with |V • B\ = to machine 
precision everywhere on the grid. 



4 ADIABATIC MOUNTAINS 

In this section, we compute Grad-Shafranov equilibria for 
several adiabatic EOS using the method described in Section 
[2] and validated in Section [3] Table [T] lists the parameters of 
each EOS, corresponding to different depth intervals within 
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Model 


K (cgs) 


r 


Equation of State 


A 


1.0 x 10 8 


1 


Isothermal 


B 


3.2 x 10 12 


5/3 


Non-relativistic degenerate electrons 


C 


4.9 x 10 14 


4/3 


Relativistic degenerate electrons 


D 


5.4 x 10 9 


5/3 


Non-relativistic degenerate neutrons 


E 


Variable 


Variable 


Piecewise polytropic 



Table 1. Numerical models of magnetic mountains with their associated 
EOS. In models A-D, the EOS is polytropic , with P(p) = Kp T ', where K 
is measured in cgs units (dyn g _r cm 3r_2 ) (1 Shapiro fc T eukolskv 1983). 
In models A-D, K and Y are held constant as M a varies. In model E, 
K and Y assume average values, which depend on M a (see Section [5]). 



the stellar crust (see Section [SJ. The scalings of the mag- 
netic dipole moment \i and mass ellipticity e with accreted 
mass M a are studied in Sections 14.11 and 14.21 respectively. 
The maximum density and local magnetic field strength are 
computed in Sections 14.31 and 14.41 respectively. In Section 
14.51 we compare the equilibrium density and magnetic field 
distributions for adiabatic and isothermal magnetic moun- 
tains. For each model in Table [TJ we stop our simulations 
once \Aip/ip\ is less than 5 per cent averaged over the grid 
(see Appendix [A]) . 

4.1 Magnetic burial: \i versus M a 

As accretion proceeds and the initial dipolar magnetic field 
lines are distorted, magnetic energy is transferred from the 
dipole to higher order multipole moments. The north-south 
antisymmetry of B r precludes the existence of even multi- 
poles. Fig. [2] displays the magnetic dipole moment \i (nor- 
malized by its initial, or surface, value) as a function of the 
accreted mass M a for models A-E in Table [TJ The maxi- 
mum accreted mass for which the iterative solver converges 
reliably (grid- averaged residual ^ 5 per cent) depends on 
the EOS, with M a , max «lx 10~ 3 , 3 x 10~ 8 , 2 x 10~ 7 , 3 x 
10 -6 Mq for models A-D, respectively. (As a corollary, the 
gradient d/i/dM a in the vicinity of the rightmost data point 
for each model in Fig. [2] is unphysically steep.) The method 
we use to calculate the dipole moment differs slightly from 
that in PM04; we integrate ip directly rather than B r , ac- 
cording to 

l(2l + l)R l m f 1 dfUgosfl) 

» i = 2(1+1) 7V ( } ^ m ' } ^(^r (14) 

for the /th multipole moment, circumventing one set of nu- 
merical derivatives and improving the accuracy of the re- 
sults. Equation (|14[) is ~ 10 per cent more accurate than 
equation (34) in PM04 for a 64 x 64 grid. The discrepancy 
shrinks to < 1 per cent for a 256 x 256 grid. 

It is clear from Fig.[2]that the characteristic mass M c re- 
quired to significantly distort the initial configuration varies 
with the EOS. If we define M c to be the accreted mass that 
halves \i from its initial value j ij, to be consistent wit h the 
empirical scaling introduced bv lShibazaki et al.1 (|l989h . viz. 

M = /Xi(l + Ma/M c )" 1 , (15) 

then Fig.[2]yieldsM c ,A « 5xlO" 4 M , M C , B « 3xlO" 8 M , 
M c ,c « 1 x 1O~ 7 M and M C , D « 2 x 1O~ 6 M for models 
A-D in Table [T] Plainly, varying the EOS makes a big dif- 
ference. M c is reduced by a factor of between 3 x 10 2 (model 



D) and 2 x 10 4 (model B) relative to an isothermal moun- 
tain. This is because adiabatic mountains are up to ~ 10 2 
times taller than isothermal ones for M a = M c (see Fig. [3] 
below and Section [475]) . At higher altitudes, the magnetic 
stress (oc r -6 ) is weaker and hence the pressure gradient 
pushes the magnetic field sideways more than in an isother- 
mal mountain. 

In the limit of small M a , one can show (see Appendix lB]) 
that the scaling of the characteristic mass M c for adiabatic 
mountains is proportional to the square of the magnetic field 
strength, as for isothermal magnetic mountains (see Section 
3.2 in PM04). Additionally, M c is also inversely proportional 
to an extra factor J(Ao,r) (evaluated as a contour plot on 
the Ao-T plane in Fig. lBl[h which depends only on the EOS 
parameters and the accreted mass through equations (|B18[) 
and (|B22[) . In this limit, one finds the following scalings of 
the magnetic dipole moment: \i oc (1 — k A M^B~ 2 ) (model 
A), fi oc (1 — kB,r>M^ 5 B~ 2 ) (models B and D) and \i oc 
(1 - k c Mi /2 B~ 2 ) (model C), where k A ,b,c,d are constants. 
We confirm in Section [5] that the realistic EOS (model E) is 
well approximated by model C and hence follows the same 
scaling. It is important to note that these /i(M a ) scalings are 
only valid in the small- M a limit (i.e. M a ^ M c , where M c is 
EOS-dependent). For M a > M c , the analytical solution no 
longer applies and numerical results have to be used. 

In Fig. [3] we plot fi/fii as a function of altitude above 
the surface for models A-D by replacing R m with r in equa- 
tion (fT4)) . The purpose is to illustrate how the screening 
currents are distributed radially for different EOS. The ac- 
creted masses are chosen to be the characteristic masses M c 
of each model in Table [1] The dipole moment turns up by 
~ 5 per cent at r ~ R m because the Neumann boundary 
condition dip/dO = 0, which holds the field lines perpen- 
dicular to the outer grid boundary, does not apply strictly 
to a dipole field. For the isothermal mountain (model A), 
the screening currents are located 10 1 — 10 2 times closer to 
the neutron star surface than in models B-D, and the iso- 
density contours contract towards the surface by the same 
factor (see Section |4~5|) . 

4.2 Mass quadrupole: e versus M a 

Fig. [4] displays the mass quadrupole moment of the moun- 
tain, expressed in terms of the mass ellipticity e, as a func- 
tion of M a . The ellipticity is given by e = \I ZZ — I yy \/Io, 
where Uj denotes the moment-of-inertia tensor, the z-axis 
lies along the magnetic axis of symmetry, and we define 
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Figure 1. Convergence between the isothermal (r = 1) and adi- 
abatic (r — > 1) magnetic mountain models at M a = 10 _5 Mq 
for the mass ellipticity e (top panel), magnetic dipole moment \i 
(middle panel) and grid-averaged residuals (bottom panel) as a 
function of adiabatic index Y . The data points represent magnetic 
mountains modelled with the EOS in equations (|11J) (red dia- 
monds), (|12|) (green rectangles) and (|13|) (blue triangles). Isother- 
mal results are denoted by filled black circles in the top two pan- 
els. 
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Figure 2. Magnetic dipole moment /i, computed at the outer 
edge of the grid and normalized to its surface value, as a function 
of accreted mass, M a , measured in solar masses, for models A 
(black diamonds), B (red circles), C (green squares) and D (blue 
triangles). Values of the characteristic masses M c are plotted as 
vertical lines for models A (solid line), B (triple-dot-dashed line), 
C (short-dashed line), D (long-dashed line), and coloured accord- 
ingly. 
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Figure 3. Magnetic dipole moment /i, calculated as a function 
of altitude (in centimetres) from the inner boundary of the com- 
putational mesh R[ n and normalized to its surface value, for the 
characteristic masses M c , measured in solar masses, for models 
A (solid black curve), B (triple-dot-dashed red curve), C (short- 
dashed green curve) and D (long-dashed blue curve). The altitude 
where the dipole moment reaches a minimum denotes where the 
screening currents end. 



Io — (2/5)M*Rf n . To zeroth order, both M a and e are pro- 
portional to the surface density p S urf- Hence, the ellipticity 
is proportional to accreted mass for M a < M c . At M a ~ M c , 
the hydrostatic pressure overwhelms the Lorentz force and 
the mountain spreads laterally, distributing the extra ac- 
creted mass evenly over a larger area (the enlarged magnetic 
polar cap) and moderating the growth of the ellipticity such 
that de/dM a < 1/M . 

The apparent turnover in e after it peaks in Fig. [4] is a 
numerical artefact, which sets in as the convergence of the 
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Figure 4. Mass ellipticity e, as a function of accreted mass, M a , 
measured in solar masses, for models A (black diamonds), B (red 
circles), C (green squares) and D (blue triangles). Values of the 
characteristic masses M c are plotted as vertical lines for models 
A (solid line), B (triple-dot-dashed line), C (short-dashed line), 
D (long-dashed line), and coloured accordingly. 

numerical algorithm worsens (see Section I4.ip . In reality, 
for M a > M c , the ellipticity saturates at the value where 
de/dM a = in a hard-surface model. IWette et all (|2010h 
examined accretion on to a non-rigid neutron star crust, 
thereby allowing the accreted matter to sink, and showed 
that the ellipticity does not saturate (i.e. de/dM a > 0) up 
to M a < 0.12Mq. Despite this monotonic increase, the el- 
lipticity of soft-surface mountains is always less than that of 
hard-surface mountains by 25-60 per cent in the mass range 

1.2 x 10" 4 < M a /M < 1.2 x 10" 1 . 

4.3 Equatorial magnetic compression 

The accreted matter transports frozen-in magnetic flux 
equatorward as it spreads sideways under its own weight. As 
a result, the magnetic field lines are 'pinched' near the sur- 
face at t he equator and flare outwards at higher altitudes lik e 
a 'tutu ' (jMelatos fe Phinnevl[200ll ; IPavne fe Melatosl 12006). 
The maximum magnetic field strength |2?| max in the equa- 
torial belt is computed as a function of M a and graphed 
in Fig. [5] for models A-D in Table [U Naturally, the lati- 
tude where \B\ maximizes moves towards the equator as 
M a increases, and the equatorial belt narrows. From Fig. [5l 
we see that adiabatic magnetic mountains produce a larger 
\B\ max (and hence a narrower belt, by flux conservation) 
than isothermal ones with the same M a . Referring to Fig. 
H this can be understood as follows. The top panel of Fig. 
[6] shows the equilibrium magnetic field configuration of an 
adiabatic and an isothermal mountain, at their characteris- 
tic masses M c (these masses are different since M c is EOS- 
dependent). At equilibrium, the hydrostatic pressure gra- 
dient at the base of the mountain (dotted red/blue arrow 
for adiabatic/isothermal mountain in Fig. [6} is balanced by 
magnetic stresses (red/blue arrow for adiabatic/isothermal 
mountain in Fig. [6} within the equatorial magnetic belt (the 
extent of the magnetic belt is denoted by red/blue shaded 
regions for adiabatic/isothermal mountains in Fig. [6}. The 




Figure 5. Maximum field strength in the equatorial magnetic 
belt at hydromagnetic equilibrium, |.B|max (in gauss), as a func- 
tion of accreted mass, M a (in solar masses), for models A (black 
diamonds), B (red circles), C (green squares) and D (blue tri- 
angles). Values of the characteristic masses M c are plotted as 
vertical lines for models A (solid line), B (triple-dot-dashed line), 
C (short-dashed line), D (long-dashed line), and coloured accord- 
ingly. Overplotted are curves of the maximum yield magnetic field 
strength, B y - ie id (in gauss), at the base of a mountain of mass M a 
for models A (solid curve), B (triple-dot-dashed curve), C (short- 
dashed curve), D (long-dashed curve), and coloured accordingly. 
Adiabatic magnetic mountains (models B— D) surpass |-B| max , de- 
forming plastically, while the isothermal mountain (model A) does 
not exceed this threshold, and hence does not crack. 

hydrostatic pressure gradients for both EOSs are compa- 
rable at characteristic accreted masses, because the mag- 
netic field lines are bent by a similar angle for all models 
at M a ~ M c [since /jl(M c ) is EOS-independent]. This can 
be expressed equivalently in terms of the comparable width 
of the equatorial magnetic belt of both mountains, since 
comparable deformation angles of the magnetic field lines 
result in corresponding widths of the magnetic belt. Refer- 
ring to the bottom panel of Fig. [6j the hydrostatic pres- 
sure gradient at the base of the accreted layer is greater 
for adiabatic mountains than isothermal ones at an equiva- 
lent M a , because M c; a > M c; d > M c; c > M c; b, where the 
subscripts A-D denote the models in Table [T] (see Section 
I4.1[) . Hence, magnetic-field lines of an adiabatic mountain 
are more deformed than those of an isothermal one to coun- 
teract this. This decreases the lateral extent of the magnetic 
belt and, by magnetic flux conservation, |I?|max increases as 
the belt shrinks. This explains why the point where |-B| ma x 
is reached moves equatorward as M a increases, and why 
|£?|max is greater for an adiabatic rather than an isother- 
mal mountain for the same M a . 

The compressed magnetic field can surpass the yield 
strength of the crust, at which point the magnetic stresses 
break the Coulomb lattice as the field deforms. Taking the 
breaking strain of the neutron star cru st to be ~ 0.1 from 
recen t molecular dynamics simulations ([Horowitz &; Kadaul 
2009), the magnetic fiel d strength at which the crustal mat- 
ter yields ([Romanilll990h is 

£ yie i d = 1.2 x 10 14 ZA- 2/3 (p/10 n g cm -3 ) 273 G, (16) 
where Z and A are the mean atomic and mass numbers, 



8 M. Priymak et al. 





Bfinal(M C ;|) 


^initial 




i5 ^^>5\ / 








NS \W>^ 




interior Vf^P*^ ^sv, \ 



Bfinal(M a; |) 




Figure 6. Schematic diagram (not to scale) of the magnetic field 
of a neutron star during magnetic burial in the case of an adia- 
batic (thick black curve) and isothermal (thin black curve) mag- 
netic mountain. An adiabatic magnetic mountain extends further 
above the surface. Undistorted dipole magnetic field lines before 
accretion are shown as dashed curves. Adiabatic and isothermal 
mountains of characteristic mass M c (top panel) and equivalent 
accreted mass M a (bottom panel) are represented [M c is the EOS- 
dependent characteristic mass defined in equation (|15|) ]. The dis- 
torted magnetic field lines during magnetic burial are shown for 
both the adiabatic and isothermal EOS of the accreted matter 
(thick red curves and thin blue curves, respectively). Subscripts 
I and A denote isothermal and adiabatic EOS, respectively. The 
Lorentz force (red/blue arrows for adiabatic/isothermal EOS) of 
the compressed magnetic field in the equatorial belt (the ex- 
tent of the belt is denoted by red/blue shaded regions for adi- 
abatic/isothermal EOS) balances the hydrostatic pressure gradi- 
ent (red/blue dotted arrow for adiabatic/isothermal EOS) at the 
base of the mountain. As more matter accretes, the increasing 
hydrostatic pressure gradient at the base of the mountain is com- 
pensated by the enlarged Lorentz force in the magnetic belt. This 
compresses the magnetic belt further towards the equator. 

respectively. We evaluate B yie id at the base of a moun- 
tain of mass M a from the n uclides present at base pr essure 
(jHaensel fe Zduniklll990aH bl: IChamel fe Haenselll2008h . The 
results are plotted as curves in Fig.[5]for the models in Table 
[1] In an isothermal mountain, we find |-B| max < ^yieid, so 
that the accreted matter does not crack and remains poly- 
crystalline, with a frozen-in magnetic field. As the substrate 
of an isothermal mountain does not spread significantly, e 
and M c are larger. Indeed, strictly speaking, crustal freez- 
ing should be included in the boundary conditions of an 
isothermal mountain calculation (implemented dynamically 
at the depth where it first occurs). On the other hand, adi- 
abatic mountains compress the magnetic field in excess of 



B y i e id for M a > 3 x 1O~ 7 M and M a > 6 x 1(T 9 M for 
models D and C respectively, while B y i e \d is surpassed for 
all accreted masses in the case of model B. This suggests 
that the accreted matt er continuously cracks or flows plas- 
tically at most depths ([Horowitz &; Kadaull2009h , validating 
the fluid approximation for models with T ^ 4/3. 



4.4 Maximum Density 

The maximum density at the base of a magnetic mountain 
is reached at the magnetic pole (see Section |43|) . We extract 
the maximum density p m ax(-Rin, 0) as a function of M a from 
the simulated models listed in Table [T] and graph the results 
in Fig. □ 

A deficiency of isothermal mountains, noted by PM04, 
is the unrealistically high density at the base, which exceeds 
the neutron drip point[j ^nd ~6xl0 11 g cm -3 at relatively 
small accreted masses of - 1O" 8 M (cf. M a - 1O _1 M 
in a typical LMXB). In contrast, Fig. [3 shows that p ma x 
is several orders of magnitude lower for an adiabatic EOS; 
none of the adiabatic mountains surpass pnd for M a < M c . 

At accreted masses approaching M a ~ 1O _2 M , models 
C and D attain crust-core densitieqj pec ~ 2 x 10 14 g cm -3 
at their bases. These models are good approximations to the 
EOS of the neutron star crust at p > 10 9 g cm -3 and p > 



10 g cm - , respectively (see Section[5]). On the other hand, 
model B does not reach the crust-core interface because it 
is too stiff and approximates the true crustal EOS only at 
low densities 10 5 < p/(g cm -3 ) < 10 7 . In the isothermal 
mountain (model A), p max exceeds pec for M a > 1O _5 M . 



4.5 Hydromagnetic structure 

A meridional cross-section of the magnetic mountain pro- 
duced by models A-D in Table [1] is displayed in Fig. [8] for 
M a = M c . The magnetic field lines and isodensity contours 
are graphed as solid and dashed curves, respectively; the 
shading also represents the density and is included to guide 
the eye. Note that the vertical scale changes dramatically 
from panel to panel. Adiabatic mountains stand 10 1 — 10 2 
times higher than an isothermal mountain for M a — M c (see 
also Section HT]) . Moreover, one finds p — »• as r — ► oo in an 
isothermal mountain, whereas an adiabatic mountain drops 
to p — at a finite altitude. Mountains with an ideal de- 
generate electron gas EOS (models B and C in Table [1]) are 
approximately 1 order of magnitude taller than those with 
an ideal degenerate non-relativistic neutron gas EOS (model 
D). 

The polar (r p ) and equatorial (r e ) mountain heights 
as well as their ratio S can be estimated analytically in 
the small- M a approximation developed in Appendix [B] For 



1 Pnd depends on wheth er the crust is cold-catalyz ed or accreted, 
as well as the exact EOS (IChamel Hae nsel 2008). Here we con- 



>ressi ble liquid drop mode l t or the hub oi the 
(lHaensel fe Zdunikl Il990al lbt IChamel fe Haensell 



sider the compressible liquid drop mode l f or the EOS of the 
accreted crust 
120081) . 

2 We choose pec large enough to contain the crust-core 
transitions of both t he Friedman-Pandh a ripande-Skyrme and 
Skyrme-Lyon EOS (|Pethick et al.l Il995l : lHaensel et all 120071 ; 
IChamel fc H aensel 2008|). 
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Figure 7. Maximum density at hydromagnetic equilibrium, p m ax 
(g cm -3 ), as a function of accreted mass, M a (in solar masses), 
for models A (black diamonds), B (red circles), C (green squares) 
and D (blue triangles). Values of the characteristic masses M c are 
plotted as vertical lines for models A (solid line), B (triple-dot- 
dashed line), C (short-dashed line), D (long-dashed line), and 
coloured accordingly. Overplotted are lines of the neutron drip 
density p^D — 6 x 10 11 g cm -3 and the density at the crust-core 
interface pec — 2 X 10 14 g cm -3 . Adiabatic magnetic mountains 
(models B-D) satisfy p m ax < Pnd> while isothermal magnetic 
mountains exceed pnd and pec f° r -&fa > 3 x 10 
1 x 1O -5 M0, respectively. 



M and M a > 



models B, C and D, and fiducial neutron star values (see 
Section [3}, we obtain 



^p|B 



5.2 x 10 4 cm, r e | B = 7.2 x 10 3 cm, S B 



7.2 x 10°, 
(17) 



r p | c = 2.8 x 10 4 cm, r e |c 



,3 x 10 3 cm, Sc = 3.3 x 10°, 
(18) 



r p | D = 6.2 x 10 3 cm, r e | D = 2.1 x 10 2 



cm 



l, Sb = 2.9 x 10 1 . 

(19) 



Comparing with Fig. [51 we see that the analytic formula in 
Appendix 151 generally overestimates r p and underestimates 
r e . This discrepancy arises because the small- M a approxima- 
tion assumes the magnetic field is nearly dipolar, whereas, 
at M c , the dipole is significantly deformed. For M a <C M c , 
there is better agreement between the numerical and ana- 
lytical solutions. 



5 CRUSTAL EQUATION OF STATE 

A realistic crustal EOS is not a simple polytrope. It in- 
cludes various pressure contributions from thermal elec- 
trons, relativistic/non-relativistic degenerate electro ns, non- 
relativistic degenerate neutrons and the ionic lattice (|Brownl 
2000). These partial pressures depend on the composition of 
the crust; accreted matter undergoes nuclear reactions (e.g. 



electron captures and pycnonuclear fusion) as the mass den- 
sity and electron F ermi energy of the compr essed matter in- 
creases with depth (|Chamel &; Haensell 20081 ). Our models of 
magnetic mountains in accreting X-ray systems necessitate 
the inclusion of a realistic accreted EOS of the neutron star 
crust. In this section, we start fro m the realistic crusta l EOS 
investigated by other autho rs (|Negele &; Vautherinl Il973l : 
|Paczvnskil[l983l : iBrownl lioooh and derive an equivalent ef- 
fective adiabatic EOS (if e ff,r e ff) as a function of M a . This 
EOS is labelled model E in Table[Q The magnetic mountains 
produced by this more realistic EOS are compared with the 
pure adiabatic ones from Section [H 

We adopt the one- component plasma appro xi mation 
for the accreted matter (["Haensel &; Zduniklll990al lbl: IBrownl 
l2000l : IChamel fe Haensell 1 20081). together wit h the n uclear 
composition proposed by lHaensel Zdunikl (|l990al lbh at 
temperature T — 10 8 K. This temperature is representative 
of the steady-state thermal profile for 10 6 < p/(g cm -3 ) < 
10 14 in accreting neutron stars containing no exotic mat- 
ter such as a pion condensate or s trange quarks in their 
interior (Mir alda-Escude et al.l Il990h . The foregoing as- 
sumptions hold for accretion rates in the range —11 < 

log lo [M/(M yr" 1 )] < -10. 

Recent w ork by I Read et ah! (|2009l ) [see also 

IVuille &: Ipserl (Jl999)] produced a four-parameter fit 
to the set of candidates for high-density EOSs in order to 
systematize the study of various observational constraints 
on the EOSs. The low-density EOS was assumed to be that 
of gro und-state cold matter given by iDouchin &; Haensel 
(2001), while the high-density candidate EOSs were 
parametrized by three free-parameter piecewise poly tropes. 
Since the EOS of th e accreted crust is stiffer than that of 
a cold-catalyzed one (|Chamel fe Haensell |2003 ). making the 
radius of a 2— IMq star 50-200 m larger than that in the 



cold-catalyzed case ([Zdunik &; HaenseTl201lh , the effective 
polytropic form for the accreted crust calculated in this 
section can be combined with observations of accreting 
neutro n stars to constrain the parameters of the parametric 
EOS of iRead et all (|2009h . 

5.1 Partial pressures 

There are three principal contributions to the pressure in 
a mountain at densities p < pec- They are as follows, 
(i) Electron pressure: this is exerted by non-r elativistic, 
relativistic and thermal electron populations ([Paczvnskil 
1983). (ii) Lattice pressure: the ionic lattice exerts nega- 
tive pressure due to electrostatic interactions within the 
Wigner-Seitz cells. It is calculated by fitting to the free 
energy in Monte Carlo simulation s of a one-component 
plasma ([Farouki &; Hamagu chi 1993). (hi) Neutron pressure: 
the effect is included in a cold-catalyzed EOS which is 
parametrized to fit 11 ground state nuc lei above the neu- 
tron drip line ([Negele &; Vautherinl \1 973). We sum the par- 
tial pressures (i)-(ii) subject to pressure continuity across 
reaction surfaces, matchin g to the cold-catalyzed EOS of 
Nege le &; Vautherinl (| 19731 ) at densities above neutron drip. 
Although the ground-state and accreted crusts contain dif- 
ferent nuclei, their respective EOS are indistinguishable for 
p > 10 13 g cm -3 , where the composition-insen sitive neutron 
pressure dominates ([Chamel &; Haen sel 2008). 

The pressure and the adiabatic index F = 
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Figure 8. Hydromagnetic structure of axisymmetric mountain equilibria at M a = M c , showing magnetic field lines (solid 
blue curves) and isodensity contours (dashed black curves) for model A with M a = 5.2 x 1O -4 M0 (top-left panel), model 
B with M a = 2.8 x IO^Mq (top-right panel), model C with M a = 1.2 x 1O _7 M (bottom-left panel) and model D 
with M a = 2.0 x 1O -6 M0 (bottom-right panel). Density contours are drawn for ?7Pmax (maximum at the pole), with 
Pmax,A = 3.9 x 10 15 g cm" 3 , p max ,B = 6.4 x 10 8 g cm" 3 , p max , c = 7.4 x 10 9 g cm" 3 , p max ,D = 3.8 x 10 11 g cm" 3 , 
r] A = 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 10~ 2 , lO" 3 , 10~ 4 , 10~ 6 , 10" 8 , lO" 10 , lO" 14 , and ry B ,c,D = 
0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 10~ 2 , 0. 



d(logP)/d(logp) of the resultant EOS are graphed 
versus density in Fig. [9] Although some parts of the EOS 
are piecewise adiabatic, other parts are not. At certain 
densities where electron capture reactions occur rapidly, 
e.g. p > 10 9 g cm -3 , the density jumps discontinuously to 
compensate for the sharp decline in electron pressure at a 
compositional interface. This behaviour is accompanied by a 
sharp drop in the adiabatic index. These discontinuities are 
an artefact of the one-component plasma approximation. 
The presence of nuclear reactions softens the EOS for 
10 12 < p/(g cm -3 ) < 10 13 , relative to uniform composition, 
whereas the addition of neutron pressure stiffens the EOS 
for p > 10 13 g cm -3 . 



5.2 Effective polytrope 

The realistic EOS [K(p),T(p)] in Fig. [9] is transformed into 
an effective adiabatic EOS, of the form P = K e fip Feff , by 



computing the mass- weighted averages 



fdr r 2 p(r)K(p) 
f dr r 2 p(r) 



fdr r 2 p{r)Y{p) 
f dr r 2 p(r) 



(20) 



for a spherically symmetric accreted layer of mass M a whose 
density profile p(r) satisfies hydrostatic equilibrium. For 
simplicity, we ignore general relativistic effects and assume 
the acceleration due to gravity to be uniform, as in models 
A-D. 

The scaling of K e ^ and r e ff with M a is shown in Fig. 
1101 The large radial variations of K and Y within the crust 
imply that K e s and r e ff depend strongly on the maximum 
achieved density and hence M a . The mass- weighted aver- 
ages are dominated by the base of the mountain. Hence, 
the mass ellipticity and magnetic dipole moment of an adia- 
batic mountain with a realistic nuclear EOS depend on M a , 
not just through the weight of the accreted layer and the 
confining magnetic stresses but also through the density- 
dependent thermodynamics at the mountain's base. 
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Figure 9. Pressure as a function of density (left-hand panel) and adiabatic index as a function of density (right-hand panel) 
for the full nuclear EOS (solid blue curve), minus neutron pressure (long-dashed blue curve) and minus compositional variations 
(short-dashed blue curve). For comparison, the isothermal EOS (triple-dot-dashed red line) is also plotted. The neutron drip and 
crust-core interface densities are marked with vertical dotted lines. Discontinuities in the adiabatic index for p > 10 9 g cm -3 
are caused by density jumps between compositional layers. The axes are log-log and linear-log in the left- and right-hand panels, 
respectively, with pressure and density measured in units of erg cm -3 and g cm -3 , respectively. 



We simulate magnetic mountains with the EOS of 
model E for 10" 10 < M a /M < 10" 7 by utilizing K eS (M a ) 
and r e ff(M a ) given in Fig. [TO] The dipole moment /x, ellip- 
ticity e, maximum magnetic field strength |-B| ma x and max- 
imum density p ma x of model E are compared with those of 
models B, C and D in Fig. [TU The hydromagnetic equilib- 
rium for model E is also plotted in Fig. [TT] for M a = M c 
(cf. Fig. [8]). As the partial pressures are dominated by rel- 
ativistic degenerate electrons at 10 -8 < M a /MQ < M c , 
fi, e, | B | ma x and p ma x in model E behave like in model C 
at these accreted masses. At M a < 10~ g Cm , -D max 
of model B approaches that of model C due to the pres- 
ence of non-relativist ic degenerate electron gas. Relativist ic 
degenerate electrons dominate model E for M a ^ M c , so 
model C can be used to approximate the realistic mountain 
for M a < 1(T 2 M . For M a > 1(T 2 M (e.g. in LMXBs), 
the dominant partial pressure comes from degenerate non- 
relativistic neutrons (model D). 

For the rang e of natal magnetic fields 1 12 < B*/(G) ^ 
10 13 inferred bv lFaucher-Giguere fe Kasoil (|2006h . M c stays 
within the range where model C applies [see equa- 
tions (lB26l) - (lB29|) of Appendix Q5J which give the model- 
dependent scaling of M c with respect to B and M a ]. If B* 
rises to ~ 10 15 G, appropriate for a magnet ar, M c increases 
from 1.2 x 10~ 7 to 2.6 x 1O~ 4 M for model C, still below 
where degenerate neutron pressure dominates (see Fig. fT0|) . 
We can therefore use model C to calculate M c under all 
plausible astrophysical scenarios. 



6 APPLICATION TO GRAVITATIONAL- WAVE 
SPIN STALLING 

Several mechanisms can brake the spin-up of an accret- 
ing neutron star: the magnetospheric centrifugal barrier 
(llllarionov Sz Sunvaevl 1975; Ghos h fe Larnbl 1 19791 ). GW 
emission (Wagoner 1984; Bildsten 1998) and the magnet ic- 
dipole torque (jOstriker &; Gunnlll969l ). Every one of these 



mechanisms eventually balances the accretion torque and 
stalls the spin-up process, when the spin frequency is s is 
large enough. We use equation (|2ip for spi n balance whic h 
assumes the usual thin-disc accretion model (Bildsten 1998). 
It should be noted that this is not necessarily valid, as 
more refined accretion models weaken the spin-up torque 
or strengthen the propeller effect, thus obviating the need 
for a strong GW torque. The feedback provided by radia- 
tion pressure in rapidly accreting systems could lead to a 
thick and sub-Keplerian inner accretion disc, which mod- 
ulates the accretion to rque of the standard thin-disc model 
(jAndersson et al1 [2005). Also, for weak accretors, if the mag- 
netospheric radius becomes larger than the corotation ra- 
dius, the s tar can exist in either a strong or weak 'propeller' 
phase fsee iRomanova et al.1 (|2008l ) and references therein), 
with the transition between these phases being strongly de- 
pendent on the kinemati c viscosity and magnet ic diffusivity 
of the accreting matter ((Romanova et al .11 20041 . 1 20051 ). Nev- 
ertheless, these improved accretion models do not invalidate 
any of the proposed GW- generating mechanisms. 

In this section, we investigate how the stalling frequency 
depends on the EOS, if all the braking comes from gravita- 
tional radiation reaction. In this work, we do not consider 
radiation-pressure feedback on the accretion disc since we 
are interested in modelling moderately accreting LMXBs 
where this effect is small. Also, in the vicinity of the bottom 
magne tic fi eld [10 7 — 10 8 G; se e Ivan den Heuvel &; Bitzarakil 
(1995) and IZhang fe Koiimal ((2006)], where the magneto- 
sphere touches the stellar surface and the propeller effect 
can be neglected, the GW torque dominates the magneto- 
centrifugal and magnetic-dipole torques. Clearly, this ap- 
proach yields an upper bound on z/ s ; the other mechanisms 
can lower v s further. 

We synthesize five Monte Carlo populations of LMXBs, 
whose spins are such that their gravitational radiation re- 
action torque exactly balances the accretion torque. We 
assume that each simulated LMXB population undergoes 
magnetic burial according to one of the five EOS in Table 
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Figure 10. Effective adiabatic EOS for mountains of different sizes (model E in Table X e ff m c g s units (left-hand panel) and 
adiabatic index T e ff (right-hand panel) as functions of accreted mass M a (in solar masses). 



\T\ The number of neutron stars in each population is chosen 
large enough (~ 10 5 ) to yield an accurate cumulative spin 
distribution. We assume fiducial neutron star parameters 
(see Section [3]) and solve 



2.09 xlO 1 Hz 



M 



1O- 1O M yr- 1 



1/5 



io- 



-2/5 



(21) 



for the equilibrium spin frequency, assuming the wob ble an- 
gle a tends to a — tt/2 due to GW back reactio n ([Cutler! 
2002) or crust-core coupling (|Alpar fe Saulisll 19881 ). The ac- 
cretion rates are selected from th e empirical luminosi ty func- 
tion of Galactic LMXB sources (jGrimm et al.| [2002). 



N(> L) = 105 



10 36 erg s" 1 



10 36 erg s- 1 



(22) 



where L is the apparent luminosity in the 2 — 10 keV 
band, and L ma x is the cut-off luminosity, combined with 
the luminosity-dependent mass fraction of the Galaxy which 
is visible to the R XTE All-Sky Monitor [see fig. 11 of 
iGrimm et al.l (|2QQ2h ] . The long-term average bolometric lu- 
minosity is related crudely to the accretion rate by the fa- 
miliar expression. 



M « LR*/GM*. 



(23) 



The results of the Monte Carlo simulations are shown 
in Fig. 1121 where we compare the cumulative distribution 
function of our spin-equilibrium models with the observed 
distribution of nuclear- powered millisecond pulsars (NMPs) 
(i.e. sources that show brightness oscillations in the tails 
of Type I X-ray bursts) , accretion- powered millisecond pul- 
sars (AMPs) (i.e. sources that exhibit X-ray pulsations) and 
accreting millisecond X-ray pulsars (AMXPs) (i.e. sources 
that exhibit either millisecond burst oscillations, X-ray pul- 
sations or both). We o btain data o n the spins of these 
objects from table 1 of IWatts etHI (|20Q8h . To be consis- 
tent wi th contemporary literature on millisecon d X-ray bi- 
naries ([Chakrabartv et al ] l2003l : lGallowavll2008h . we adopt 
the following naming convention for these sources: accret- 
ing millisecond pulsars are AMPs, burst oscillation sources 



are NMPs, and we combine these two populations into 
AMXPfl To distinguish between the confirmed and un- 
confirmed sources, we plot all/confirmed NMPs (thin/thick 
triple-dot-dashed green lines), AMPs (thick orange line) 
and all/confirmed AMXPs (thin/thick dashed blue lines). 
Curves represent cumulative distribution functions of mod- 
els A (dot-dashed black curves), B (triple-dot-dashed red 
curves), C (short-dashed green curves), D (long-dashed blue 
curves) and E (solid purple curves). We update the spin o f 
EXO 0748-676 from 45 to 552 Hz ([Galloway et al.ll2010l ). 
and we do not discriminate between intermittent pulsars 
and AMPs (i.e. those sources which exhibit intermittent or 
persistent X-ray pulsations during outburst, respectively). 

The luminosity function is denned for the RXTE All- 
Sky Monitor catalogue (2 — 10 keV band), w hich is flux- 
limited below ~ 10 35 erg s" 1 ([Grimm et al.ll2002l ). Two max- 
imum luminosity cut-offs are investigated, namely L max = 
2.7 x 10 38 erg s _1 (to include the most luminous LMXB Sco 
X-l) and 3.2 x 10 37 erg s" 1 (most luminous AMP Aql X- 
1), encompassing the luminosity range of all confirmed and 
unconfirmed AMXPs. All sources are assumed to follow the 
same power-law scaling of the luminosity function. 

The All-Sky Monitor underestimates the true bolomet- 
ric luminosity, and hence the accretion rate, due to the pres- 
ence of signi ficant hard X -ray tails > 10 keV in LMXB X- 
ray spectra ([B arret 2001). Although this can be corrected 
(Gal loway et al.l 12008). we do not attempt to do so here, 
because equation ([23]) is approximate anyway, equation ([2T]) 
depends weakly on M, and the bolometric correction factors 
differ by up to ~ 40 per cent between sources. 

Considering typical LMXB lifetimes of ~ 10 8 yr 
( Podsiadlowski et al. 2002), the accreted masses in these 
systems are evaluated to be in the range of 10 -4 < 
M a /MQ < 10 1 . Therefore, enough matter has been trans- 
ferred in these systems to reach the characteristic masses 
and saturation ellipticities for the models in Table [1] given 



3 By contrast, the naming convention used bv I Watts et al.l (2008) 
reflects how the spins are inferred observationally: accreting mil- 
lisecond pulsars, burst oscillation sources and quasi-periodic os- 
cillation sources. 
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Figure 11. Magnetic dipole moment fi, calculated as a function of altitude (in centimetres) from the inner boundary of the compu- 
tational mesh R[ n and normalized to its surface value (top-left panel), for the characteristic masses M c , measured in solar masses, for 
models B (triple-dot-dashed red curve), C (short-dashed green curve), D (long-dashed blue curve) and E (solid purple curve). Also 
plotted are: normalized magnetic dipole moment \i (top-right panel), mass ellipticity e (middle-left panel), maximum field strength 
| B | max in gauss (middle-right panel) and maximum density p m ax in units of g cm -3 (bottom-left panel), as a function of accreted 
mass M a , measured in solar masses, for models B (red circles), C (green squares), D (blue triangles) and E (purple crosses). Values 
of the characteristic masses M c are plotted as vertical lines for models B (triple-dot-dashed line), C (short-dashed line), D (long- 
dashed line), E (solid line), and coloured accordingly. Overplotted in the bottom-left panel is the line of the neutron drip density 
Pnd = 6 x 10 11 g cm -3 . The hydromagnetic structure of model E is displayed in the bottom-right panel at M a = M c = 1.5 x 1O _7 M0, 
showing magnetic field lines (solid blue curves) and isodensity contours (dashed black curves). Density contours are drawn for 77p max 
(maximum at the pole), with p max = 1.0 x 10 10 g cm" 3 , and rj = 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 10~ 2 , 0. 
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initial magnetic fields of 10 12-5 G. Hence, for each simulated 
LMXB population, we assign the ellipticities of the neutron 
stars to be the saturation values for the respective EOS in 
Tabled] 

From Fig. [12] we see that an isothermal mag- 
netic mountain (model A) stalls the star at z/ s ~ 
1 Hz (B*/10 12 - 5 G)~ 4/5 , where the B* scaling follows from 
M c oc B 2 of equation (30) in PM04 and equation fl3TJ. One 
would therefore need B* « 10 10 G to fit the observed spin 
distribution, con tradicting popu l ation synthesis studies of 
isolated pulsars jH artman et al. 1993; lArzoumanian et al.l 
2002: [Faucher-Giguere &; K asui 2006j). Adiabatic magnetic 
mountains (models B-E) are generally in better agreement 
with the observed spin distribution. In fact, models B, C 
and E produce a good fit to all of the observed spin distri- 
butions. Equation (|B26|) in Appendix [Bl for M C (B*) implies 
v s oc B* 4 ^ 9 for models B and D, and v 3 oc B* 8 ^ 15 for model 
C. Thus, a better fit to the empirical spin distributions can 
be obtained for models B and C if the fiducial magnetic 
field in the range of 10 12 - 10 13 G, rather than 10 12 5 G, is 
considered. Although model D cannot match the observed 
spin distribution in this range, it is possible that Ohmic 
diffusion can improve the agreement by allowing the moun- 
tain to spread, resulting in a lower saturation ellipticity and 
hence higher equilibrium spin frequencies. 

It appears that the equilibrium spin frequencies of con- 
firmed NMPs are systematically higher than those of AMPs; 
their cumulative distributions are offset to the right and left 
of the AMXP distribution, respectively (see Fig. [12]). This 
is qualitatively consistent with the GW spin stalling mech- 
anism, as the median time-averaged accretion luminosities 
of NMPs are ~ 20 times higher than those of AMPs, re- 
sulting in higher equilibrium spin frequencies by a factor 
of w (20) 1//5 « 1.8 (under the assumption of similar ellip- 
ticities in these systems). This roughly corresponds to the 
frequency separation between the observed NMP and AMP 
distributions in Fig. [12] supporting the GW spin stalling hy- 
pothesis. On the other hand, if outburst luminosities of these 
objects are considered instead, the separation in predicted 
equilibrium spin frequencies becomes negligible. 

Another noteworthy feature of Fig. [12] is the steep gra- 
dient of the observed distribution at v s ~ 500 Hz (D. K. 
Galloway, private communication). The theoretical curves 
for models A-E can reproduce the shape of the distribution 
for is s < 400 Hz. For the range of L max investigated, theo- 
retical curves do not rise steeply enough to fit the higher- 
frequency (y s > 400 Hz) end of the distribution. This is a 
problem for stalling models in general, not just magnetic 
mountains; the e _2//5 scaling in equation (|21|) is too gentle. 
The observed steepening could be caused by differences be- 
tween the luminosity functions of Galactic LMXB sources 
and AMPs/NMPs. Allowing for a realistic distribution of 
saturation ellipticities (e.g. due to a lognormal natal mag- 
netic field distribution of isolated pulsars, predicted by pop- 
ulation synthesis studies) worsens the steepening problem, 
if the luminosity function is assumed to be independent of 
the magnetic field. It is possible that another mechanism 
(such as the 'propeller' effect) sets z/ s , but its dependence 
on underlying variables (i.e. v s oc B* 6 ^ 7 M 3 ^ 7 ) is even gen- 
tler than gravitational radiation reaction. We defer a full 
investigation of this puzzle to a future paper. 



7 DISCUSSION 

Magnetic burial in accreting neutron stars has several im- 
portant astrophysical consequences. It creates a significant 
mass quadrupole moment, which potentially stalls the spin- 
up of an LMXB by gravitational radiation reaction. It also 
reduces the magnetic dipole moment, in accord with the 
observed \i versus M a relation in neutron st a r bin aries 
presented in fig. 2 in Taam fc van den Heuvell |l986). In 
the context of the statistical evidence agai nst field d ecay 



over 10 — 10 yr in isol ated pulsars ( Bhattach arva et al.l 
1993; lLorimer et aHll997l \ magnetic burial can be invoked 
to explain both the low magnetic fields in LMXBs and 
millisecond pulsars (|Chanmugaml Il992l ; lLamb fe Yu|[2005; 
IZhang &; Koiima| 2006) and the obs erved spin distribution 
of LMXBs (jChakrabartv et alJl2003h . However, before mag- 
netic burial is deemed a viable explanation for the above 
phenomena, the effect of the EOS on the burial process must 
be quantified. 

In this paper, we show that the effect of the EOS is 
large. Magnetic burial is more effective for 4/3 ^ V ^ 5/3 
than for r = 1, in the sense that less matter must be accreted 
in the former case than in the latter in order to achieve the 
same amount of magnetic dipole screening. For the EOS 
listed in Table Q] M c decreases from 5.2 x 1O _4 M (model 
A) to 2.8 x 1O" 8 M (model B), 1.2 x 1O" 7 M (model C), 
2.0 x 1O" 6 M (model D) and 1.5 x 1O" 7 M (model E), 
for B — 10 12 ' 5 G. Likewise, the saturation ellipticities de- 
crease from 3.2 x 10" 4 (model A) to 1.3 x 10" 8 (model B), 
6.0 x 10~ 8 (model C), 9.0 x 10~ 7 (model D) and 7.3 x 10~ 8 
(model E). This is a general result, applicable to a vari- 
ety of scenarios where magnetic confineme nt of accreted 
matt e r can occur, such as T Tauri stars (Ber tout et al.l 
1988; lHartmann et al.l Il998h, young neutron stars accret- 
ing fr om a fallback disc ([Chatteriee et al. 2000l; I Wang et all 



l2Q06h and magnetic white dwarfs (|King fe Lasotaril97fl : 
Wic kramasinghe &; Ferrariol l2000h . The characteristic mass 
scales quadratically with the magnetic field strength in all 
models but with different powers of the accreted mass: we 
have M c oc B 2 M? , with (3 = (model A), -4/5 (model B), 
-1/2 (model C) and -4/5 (model D) (see Appendix© . The 
maximum density at the base of an adiabatic mountain sat- 
isfies /?max <C 10 14 g cm -3 , unlike for isothermal mountains, 
where it is unrealistically high. We find that crust al cracking 
occurs as burial proceeds, because the yield magnetic field 
strength is typically surpassed in non-isothermal models. 

A Monte Carlo analysis of neutron stars in LMXBs, 
with M drawn from an empirical distribution and B* set 
to the fiducial 10 12 ' 5 G, shows that models B, C and E 
yield 100 < i/ s /(Hz) < 600 within the gravitational spin- 
equilibrium scenario (Bildsten 1998). This is in accord with 
the w 180-620 Hz confirmed spins of AMXPs. Model D 
predicts 50 < v s / (Hz) < 300, slightly too low to explain 
the data. In comparison, the isothermal magnetic mountain 
(model A) does not agree with the data at all, yielding is s 
values 1 order of magnitude lower than those from model D. 

We compute the magnitude of the GW strain ho of 
the AMXPs and quasi-periodic oscillation (QPO) sources 
by apply i ng the gravitational spin-equilibrium argument of 
iBildstenl (1998) to the sources in table 1 in IWatts et al.l 
(2008). Here, we differentiate between the confirmed and 
unconfirmed sources, as well as AMPs, NMPs and sources 
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Figure 12. Observed cumulative spin distributions of all NMPs (thin triple-dot-dashed green line), confirmed NMPs (thick triple- 
dot-dashed green line), AMPs (thick orange line), all AMXPs (thin dashed blue line) and confirmed AMXPs (thick dashed blue line), 
together with the theoretical distributions predicted by model A (dot-dashed black curves), model B (triple-dot-dashed red curves), 
model C (short-dashed green curves), model D (long-dashed blue curves) and model E (solid purple curves). The theoretical curves are 
based on the GW stalling mechanism of Bildsten (1998), with the saturation ellipticities e(M a = M c ) for each of the models in Table 
1. Two luminosity cut-offs are considered: L max = 2.7 x 10 38 erg s _1 (rightmost theoretical curves) and L max = 3.2 x 10 37 erg s _1 
(leftmost theoretical curves). The right-hand panel zooms into the range 90 ^ i/ s /Hz ^ 1200 in the left-hand panel and displays the 
unconfirmed sources as well. 



that exhibit both persistent pulsations and burst oscilla- 
tions. The results for AMPs (orange diamonds), confirmed 
NMPs (teal squares), unconfirmed NMPs (unfilled squares), 
QPOs (yellow triangles) and sources exhibiting both pulsa- 
tions and burst oscillations (teal diamonds) are shown on a 
wave strain ho versus wave frequency / plot in Fig. [13] where 
/ = 2z/ s . The highest / value considered here corresponds to 
2^ s ,max, where i/ S)max = 760 Hz is the maximum inferred spi n 
in NMPs via Bayesian analysis ([Chakrabartv et alj|20pj ). 
When computing ho, we assume that the transient sources 
are in torque balance during outburst. This is in accord with 
lHartman etaD (|2008T ). who argued that SAX J1808.4-3658 
is secularly spinning down between outbursts and is thus 
likely to be in spin equilibrium during outburslQ. 

The charac te ristic GW strain ho [denned in 
Ijaranowski et all (1998)] detectable by Laser Interfer- 
ometer Gravitational Wave Observatory (LIGO) and 
the proposed Einstein Telescope from a periodic source 
at a distance o f 3 k pc [representative of Sco X-l; see 
iBradshaw et al.l (l999)] with a false alarm rate of 1 per 
cent and a false dismissal rate of 10 per cent for a compu- 
tationally feasible integration time of 14 days is overplotted 
in Fig. El for LIGO S5 (thin solid curve), LIGO S6 (thin 
short-dashed curve), Advanced LIGO in the broad-band 
configuration (thin dot-dashed curve), lower envelope of 
Advanced LIGO in the narrow-band configuration (thin 
triple-dot-dashed curve) and the proposed co nventionale 
Einstein Telescope (thin long-dashed curve) (jHild et al.l 



4 The transition between the spin-up and spin-down episode 
within the 2002 o utburst of SAX J1808. 4-3658 found by 



l201ll : IWatts et al.1 120081 : ISmith et al.1 l2009h . We also plot 
ho versus / for neutron stars with magnetic mountains 
at a distance of 5 kpc, with magnetic field of 10 12-5 G, 
for models A (thick dot-dashed black curve), B (thick 
triple-dot-dashed red curve), C (thick short-dashed green 
curve), D (thick long-dashed blue curve) and E (thick solid 
purple curve). 

Model A significantly overestimates ho with respect to 
both the inter ferometer sensitivity curves and the inferred 



Bildsten (1998) limits. In contrast, model E undercuts the 



Burderi et al . (2006) is probably due to pulse shape change 
5 The xylophone configuration of the Einstein Telescope closely 
matches the sensitivity o f the conv entional configuration at fre- 
quencies > 30 Hz dHild et alJl2010h . 



Bildsten! (|1998f ) limit for QPO sources, implying either the 
natal magnetic fields of these sources are ~ 10 13-5 G, or that 
these objects are not in GW spin equilibrium. All the con- 
firmed AMXPs and most of the unconfirmed AMXPs are 
consistent with model E. They lie below the model E curve 
either because they have B* < 10 12 ' 5 G or because Ohmic 
diffusion prevents the ellipticity from saturating. We note 
that the current magnetic mountain models are still prelim- 
inary. Effects that have not yet been modelled faithfully in 
the context of magnetic burial may modify the saturation el- 
lipticities. Therefore, it is still premature to quantify the ab- 
solute detectability of magnetic mountains as GW sources. 

There have been two directed searches for GWs fr om 
the accreting neutron star Sco X-l ( Abbott et aDl2Q07al lbh. 
which is expected to be the strong est emitter of its class 
in the GW spin stalling scenario (Bildsten! 1 19981 ). The first, 
coherent search computed the F-statistic on 6 h of LIGO 
S2 data, coincident between the Hanford and Livingston 
interferometers. Assuming a non-eccentric orbit, it placed 
a 95 per cent confidence upper limit on the GW strain 
from Sco X-l of ho = 1.7 x 10" 22 in the 464-484 Hz fre- 
quency band, and hp 2.2 x 10~ 22 in the 604-624 Hz 
frequency band ([Abbott et al . 2007a), which corresponds 
to an upper limit on the ellipticity of the neutron star of 
e ~ 4 x 10 -4 . The second, semicoherent search performed 
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a radiometer analysis of 20 days of triple-coincidence LIGO 
S4 data. It yielded a 90 per ce nt confidence upper limit of 
4ms ~ 3.4 x 10" 24 (//200 Hz) (jAbbott et al.ll2007bh . As re- 
quired by the non-detect ion of gravitational e mission from 
accreting neutron stars (| Abbott et al.1 l2QQ7al lbh. adiabatic 
EOS reduce the GW detect ability of magnetic mountains 
below the current detection threshold of ho « 10 -23 . In 
comparison, the saturation ellipticities of ideal isothermal 
magnetic mountains of model A are above this threshold 
and should have already been detected. 

The models in this paper are not the final word on mag- 
netically confined mountains. The range of accreted masses 
investigated here is w ell below M a ^ lO^ M^, the typical 
value for an LMXB (|Burderi et alJ [1999). due to numeri- 
cal breakdown. If e truly saturates for M a ^> M c , then this 
failing is less serious for the GW applications than for under- 
standing /i(M a ), but it should be noted that the saturation 
hypothesis has not been tested rigorously for M a > 10M C 
(|Pavne fe Melatod 12004 IVigelius fe Melatodl2009aT ). A pre- 
cise calculation of mountain equilibria for an exact, depth- 
dependent nuclear EOS cannot be carried out within our 
Grad-Shafranov formulation, although a relativistic degen- 
erate electron EOS (model C) is a fair approximation for 
M a « M c . The models in this paper are constructed on 
an impenetrable and EOS- and M a -dependent surface R in 
within the crust, w hich prevents sinking past this boundary. 
IWette et al.l (|2010h showed that, for isothermal mountains, 
sinking reduces e by up to 60 per cent. In the presence of 
Ohmic diffusion, a balance is achieved after a mass Md is 
accreted (Md depends on magnetic field, temperature, ac- 
cretion rate and EOS), in which th e rate of cross-field mas s 
transport equals the accretion rate ( Mela tos &; Pavnell2005h . 
As our model is not time-dependent, the Hall effect is also 
missing. Hall d rift acts to break down the m agnetic field to 
shorter scales ([Hollerbach &; Rudigerll2002U2004l) and may 
operate in isolated neutron stars ([Rheinhardt &; Geppertl 
2002; Rhe inhardt et al.1 2004) but is thought to be relatively 
unimportant in accreting neutron stars, where it is domi- 
nated by Ohmic diffusion ([Cumming et al.ll2004h . The crys- 
talline lattice of the crust is thought to melt in thin layers 
where electro n captures h ave significantly reduced the nu- 
clear charge ([Brownl l2000h . This is expected to have non- 
negligible effects on magnetic burial, as the boundary con- 
dition on the magnetic field becomes a function of density 
rather than radius (line- tying where solid, free where liquid). 
Finally, the three-dim ensional stability of MHD equilibria 
depends on the EOS (|Kosihski fe Hanaszl l2006h . We leave 
the investigation of these phenomena to future work. 
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(vi) The resultant F* (ip) is under-relaxed with the initial 
input value F (0) (VO via = AF (0) (VO + (1 - A)F*(^). 
The under-relaxation parameter is usually set to A = 0.9. 

(vii) Fit F{^) with a \ 2 polynomial fit of the form 
X^=o ^V^ 5 where the degree of the polynomial is typically 
N p = 8. F'(ip) is computed analytically from the coefficients 
of the polynomial. 

(viii) Values of F(ip) and F'(ip) are linearly interpolated 
across the mesh and fed into the source terms in equations 
Q, ©, © and © as required. 

(ix) The Grad-Shafranov equation is solved for the in- 
termediate quantity -0* (r, 6) by an iterative Poisson solver 
that utiliz es successive over- rela xation with Chebyshev ac- 
celeration (|Pavne fe Melatos 2004). The Poisson solver halts 
when the fractional residuals are < 10 -2 everywhere across 
the grid. 

(x) The solution for ip( n \r,6) is subsequently under- 
relaxed via V (n) = B^ (n_1) + (1 - 0)V>*. The under- 
relaxation parameter is typically given by G = 0.999 for 
r = 1 for T = 1 and 9 = 0.9995 for IV 1. 

(xi) Repeat steps (v) 



(x) until convergence is achieved. 



APPENDIX A: NUMERICAL ALGORITHM 

The iterative solver described in appendix B of PM04 must 
be modified to handle an adiabatic EOS. Equations ((5} and 
J9} show that two quantities, ^(r, 9) and F(ip), must be 
solved for at each iteration step. 

F(ip) appears implicitly in both the left- and right-hand 
sides of equation J9} for r / 1, whereas it is derived explic- 
itly from ip in one go via equation (|SJ) for V = 1. Below we 
outline briefly the major steps necessary to calculate hydro- 
magnetic equilibria in the general case. 

(i) A two-dimensional mesh is defined, of size Nx x 
Ny (typically 256 x 256). The X-axis scales proportion- 
ally to logr to increase the resolution close to the sur- 
face, where most of the screening currents lie; we write 
X = \og[x exp(Lcc) + 1], with x = (r — Ri n )/xo. The F-axis 
scales as cos#, with Y = cos 6. 

(ii) The flux function i/j(°\X,Y) is initialized across the 
mesh, with ifj( \r,6) = ip*Ri n sin 2 0/r (i.e. dipole). 

(iii) Ny — 1 contours are laid down (the number is chosen 
to minimize the occurrence of grid crossings), with contour 
values spaced linearly in Y (i.e. linear in cos#). 

(iv) Equation (|10p is used to compute dM/dip along each 
contour. 

(v) For r = 1, compute F(ip) via equation ©. For 
T / 1, iteratively solve equation J9} using the previous 
iterate F^ -1 ^ (ip) as a starting point [with F^(ip) uni- 
form]. The value of F(°\ip) is chosen to guarantee that the 
term in square brackets within the integral in ©, which 
equals the density, vanishes at the edge of the computa- 
tional grid after every iteration, so that the integral in J9} is 
always well defined. The iteration is halted when the grid- 
averaged fractional residual drops below a threshold, usually 
([F (n+1) (^)-F (n) (^)]/F (n) W) < 1-Ox 10~ 6 . Convergence 
is usually achieved after ~ 20 iterations. 



APPENDIX B: ANALYTIC APPROXIMATION 
FOR THE CHARACTERISTIC MASS M c 

An analytic formula can be obtained for the characteristic 
mass M c by calculating /x(M a ) analytically in the small- M a 
limit and looking for the value of M a where \i drops to half 
its unperturbed value. 

To calculate /i(M a ) in the small- M a limit, we follow ap- 
pendix A3 in PM04 and proceed in three steps. First, we pick 
a simple form of F(ip) which linearizes the Grad-Shafranov 
equation while approximating the exact numerical result: 



(Bl) 



Secondly, we evaluate the right-hand side of the Grad- 
Shafranov equation assuming that the flux function is ap- 
proximately dipolar: 



^*#in(l -M 2 ) 



(B2) 



This is justified because the magnetic field is weakly dis- 
torted in the small- M a limit. Thirdly, we solve the Grad- 
Shafranov equation with the above source term to obtain 
the leading-order correction to ip. 

We begin by re-expressing the radial coordinate in 
terms of the fractional altitude x 



(B3) 



In a typical mountain, with height < 10 5 cm (see Section 
I4.5p . one always has x <^ 1 within the mountain. With 
equations (|B1[) - (|B3[) , the Grad-Shafranov equation ([5]) in 
the small- M a approximation becomes 



Qo 



r^ 1 /r{Q ^[i_(i_ M 2 )(1 . 



i/(r-i) 



- x )]}(r-i)/r 

(B4) 



The Lorentz force vanishes when the right-hand side of equa- 
tion (|B4[) is zero. Therefore, the maximum height of the 
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magnetic mountain as a function of latitude can be written 

as 



2\(r-i)/r 



(r - i)0o 

for fj, 2 3> x(l — fj, 2 ), and 
_ KT r [Qo4> t ] r ~ 1 

Xmax ~ [(r-i)M r 



(B5) 



(B6) 



for fi <C x/(l — x) (i.e. near the magnetic equator). (It is 
easy to check that one has x m < 1 a posteriori for typical 
parameters.) Therefore, for adiabatic magnetic mountains, 
the ratio of polar to equatorial heights is 



[QoV>* 



-(r-ir/r 



fo(r-i) 



(B7) 



Equation (|B4[) can be solved by the method of Green's 
functions. From Section 3.1 in PM04, we write 



/l POO 
dp / dr r 2 gi +1 {r,r)C\ 
-l ^i? in 



3/2 (//)Q(rV), 



21+1 



Q(r ,n) = Qo(l-M )»• S 1- 



2(j + !)(* + 2) 
(2/ + 3) ' 



(B9) 



(BIO) 



i/(r-i) 



(Bll) 



(B12) 



with r< = min(r, r') and r> = max(r, r'). The symbol 
Cf^ 2 (/x) denotes the /th Gegenbauer polynomial. The first 
few are listed for reference: C^ 2 (n) — 1, Cf^ 2 (/i) = 3/i, 
C^J 2 (\i) — (3/2) (5/x 2 — 1). Since we are interested in how the 
dipole moment is screened at large r, we assume r > r^ax al- 
ways, where r^ax is the top of the mountain. This simplifies 
the radial Green's function to 



9i(r, r) 



r '(*-i) 
(21 + l)r l 



Rin 



21 + 1 
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Away from the magnetic equator, i.e. // 2 > x' (1 — // 2 ), 
equations (|Bl)-(|B3]), (|B5j and (lB8l)- (IBlll) combine to give 



(B14) 
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f dp' [ 

Jx' 1 / 2 JO 
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with 
Ao = 



Mr- 



dx' 
A x' 



/2(r-i)/r 



V(r-i) 



(B16) 



Our goal is to calculate the dipole moment as a function 
of M a given (iBHl) and (iBTSl) . In the limit r oo, the Z ^ 1 



contributions to /i vanish, and equations (|B14[) and (|B15|) 
reduce to 



with 



3QoRf r 
2^ 
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(B18) 



x'(l 
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Contours of /(Ao,r) are plotted in Fig. IB II for reference. 

To express Qo in terms of the other variables, we sub- 
stitute equations (|10[) . (|B1|) . ()B2|) and ()B16|) into equation 
dU) to give 



[Qo(V*-V0] 



i/r 



K 1 ^ M a ^R in exp(-^ a 
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where b — ?p*/?p a is a constant that parametrizes the lat- 
eral extent of the accretion column. Equation (|B19[) is not 
strictly an equality; the linear ansatz F(ip) = Qo(ip* — VO 
is not an exact solution in the small- M a limit (see fig. 6 in 
PM04, which presents a numerical comparison). Hence, to 
evaluate Qo approximately, it is enough to integrate equa- 
tion (|B19[) through the centre of the mountain, where most 
of the mountain mass resides (i.e. along the polar flux line 
ip = 0). This has the added advantage that the resultant 
contour integral has no dependence (ds = dr, ip = 0). 
Changing variables according to equation (|B3[) , substitut- 
ing equation (|B5[) for the upper integration limit in x , and 
taking (1 + x) ~ 1 inside the integral, we arrive at 



[Qo^* 



i/r 



K l/T 
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The expression for Qo is therefore 
M a M*Gb 



(B20) 
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and hence equation (|B16|) becomes 
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(r - 1)GM, j bexp(b)GM a M* 
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(B22) 



rK 1 ^^ I 27r[exp(6) - l]R? n 
For fiducial neutron star parameters (e.g. Section [3} and the 



adiabatic equations of state B-D in Table [T] equation (|B22|) 
reduces to 
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x 10" 



M a 
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Upon substituting equations (|B17|) and (|B21|) into equa- 
tion (fT4)) (with 1 = 1, and r instead of R m ) and comparing 
with the phe nomenological burial la w fi = in(l — M a /M c ) 
postulated bv lShibazaki et al.l ((1989) in the small- M a limit, 
we obtain 

4tt[1 -exp(-fr)]V? 2 



3M*G&I(Ao,r) 
2.8 x 1Q~ 9 [1 -exp(-b)] 



(B26) 



5 



10 12 - 5 G 



i(A ,r) 



io- 



M . 



The integral in equation (|B18[) is computed for models 
B, C and D with fiducial neutron star parameters and plot- 
ted as a function of M a in Fig. IB2I For the case Ao > 1, 



equation (|B18|) reduces to the following expressions 
J B (M) = A B M 7/S + B B M 5/S + C B M 4/5 , (B27) 
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A C M 7/6 + B c M 0/0 + C C M 
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(B28) 
(B29) 



with M M a /M , A B 3.5 x 10 7 , B B -1.7 x 10 5 , C B 
57, A c = 6.9xlO- 3 ,Bc = 2.2xl(T 2 ,Cc = 2.7xl0~ 2 , A D = 
7.2 x 10" 3 , B B = 2.0 x 10" 2 and C D = 2.7 x 10" 2 . 



In contrast to equation (|B26[) . the scaling of M c for 
isothermal magnetic mountains [from equations (29) and 
(30) in PM04] is 
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Figure Bl. Contour plot of the integral J(Ao,r) defined by 
equation (|B18|) . with Ao denned by equation (|B16|) . T and Ao 
span the typical range expected in magnetic mountains. The 
contours correspond to (from top to bottom) log 10 /(Ao,r) = 
-10, -9, -8, -7, -6, -5, -4, -3, -2, -1.7, -1.6. 
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Figure B2. Plot of equation (B18) for models B (solid red curve), 
C (dashed green curve) and D (dot-dashed blue curve), as a 
function of accreted mass M a (measured in solar masses). In the 
regime where the maximum height of the magnetic mountain is 
x < 1 (i.e. Ao > 1), /(M a ) is polynomial. I (Ma) saturates at 
« 10 -1 - 5 in the case where Ao < 1, which is an artefact of the 
approximations used in the small- M a limit. 



